

% function takes image of vessel centers selected in FIJI and draws
% lines for linescans from it; input is just either the filename of the
% manually selected image or the actual image; outputs are ep1 and ep2
% which are variables where the first column is the x-coordinate and second
% column is the y-coordinate of the first endpoint (ep1) and second
% endpoint (ep2) of the lienscan line

function [ep1,ep2,endpointIDs] = linescansFromCenters(sourceImage, preselectedCenters)

if nargin > 1 && ~isempty(preselectedCenters)
    % if already have centers selected (i.e. through automatic detection)
    % then use that
    im = preselectedCenters;
elseif ischar(sourceImage)
    % load vessel selections from FIJI:
    im = double(imread(sourceImage));
else
    im = sourceImage;
end
vesselIDs = im;          %this will keep vessel IDs if they're available
%im(im<max(im(:))) = 0;
im(im~=0) = 1;

ind = find(im == 1);
orientIm = nan(size(im));
for n = 1:length(ind)
    temp = getLocalOrientation(im,ind(n));
    if ~isnan(temp)
        orientIm(ind(n)) = temp;
    end
end
orientIm = deg2rad(orientIm);
%% downsample the number of points
downsampleDistance = 10;
[~,downsampledInd] = downsampleConnected(im,downsampleDistance);

%% draw linescans for all these points:
maxradius = 30;     %in pixels
ep1 = zeros(length(downsampledInd),2);
ep2 = ep1;
[M,N] = ind2sub(size(im),downsampledInd);
endpointIDs = zeros(size(ep1,1),1);             %saves the vessel# from which that endpoint came
for n = 1:length(downsampledInd)
    ep1(n,1) = sin(orientIm(downsampledInd(n)))*maxradius + N(n);
    ep1(n,2) = cos(orientIm(downsampledInd(n)))*maxradius + M(n);
    
    ep2(n,1) = -sin(orientIm(downsampledInd(n)))*maxradius + N(n);
    ep2(n,2) = -cos(orientIm(downsampledInd(n)))*maxradius + M(n);
    endpointIDs(n) = vesselIDs(downsampledInd(n));
end

%% since we're excluing center points that are too close to the edges, these will be NaNs in the results, so we should remove them:
ind = ~isnan(ep1(:,1));
ep1 = ep1(ind,:); ep2 = ep2(ind,:);

%% testing block to check that linescans are drawn properly
if 1
    figure
    imagesc(orientIm)
    hold on
    for n = 1:size(ep1,1)
        plot([ep1(n,1),ep2(n,1)],[ep1(n,2),ep2(n,2)],'o-')
    end
    axis equal
end



